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Abstract 

The technique of obtaining second order, oscillation free, total variation diminishing 
(TVD), scalar difference schemes by adding a limited diffusive flux (“smoothing”) to a sec- 
ond order centered scheme is explored. It is shown that such schemes do not always converge 
to the correct physical answer. The approach presented here is to construct schemes that 
numerically satisfy the second law of thermodynamics on a cell by cell basis. Such schemes 
can only converge to the correct physical solution and in some cases can be shown to be 
TVD. An explicit scheme with this property and second order spatial accuracy was found 
to have an extremely restrictive time step limitation (At < Ax 2 ). Switching to an implicit 
scheme removed the time step limitation. 

1. Introduction 

The numerical simulation of fluid flow is an interesting problem. For many cases of in- 
terest the equations and constitutive relations are well known and have been for some time. 
A great deal of excellent work has been devoted to solving these equations. The current 
state of the art for convection dominated problems is represented by various total varia- 
tion diminishing (TVD) schemes and the closely related essentially non-oscillatory (ENO) 
schemes. These schemes can be constructed using geometric arguments. Unfortunately, it is 
possible to construct wiggle free schemes that give physically impossible answers. Another 
approach, presented here, is to construct schemes that numerically satisfy the second law 
of thermodynamics on a cell by cell basis. Such schemes do not allow unphysical solutions 
and in some special cases can be shown to be TVD. 

Any general numerical method for the simulation of fluid flow must deal with the 
case of convection dominated flows with shock waves. Such cases are highly nonlinear and 
represent a real challenge to many algorithms. This work deals with such problems. 

Recently Sweby [l] proposed a method for constructing TVD schemes. It was the 
investigation of this method that led to the work of this paper, which shows how to construct 
entropy satisfying, second order accurate, central difference schemes stabilized by a limited 
diffusive flux. In some cases these can be shown to be TVD. 

All the schemes reviewed by Sweby (Roe[2],Van Leer[3], Osher-Chakravarthy[4]) were 
based on the second order Lax-Wendroff scheme [5] and consequently had the usual draw- 
backs of such a scheme. Specifically, such schemes have a steady state solution that depends 
on the time step and are, therefore, unsuitable for implicit methods. Also, not all TVD 
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schemes converge to the correct solution. It is possible to construct, for the wave equation, a 
TVD scheme that violates entropy, thus resulting in a square wave for arbitrary initial data. 
The objective of this work is to construct a second order accurate space difference scheme 
that satisfies a cell entropy inequality and converges to a solution that is independent of 
the time step. 

For entropy satisfying schemes using explicit Euler time advance, extremely restrictive 
time step limitations are needed to achieve second order accuracy. On the other hand, using 
implicit Euler time advance, a scheme can be constructed that has all the desired properties 
and no time step limitation. 

The approach presented here is still unproven but shows great promise. It is one 
dimensional though the extension to multidimensions is straightforward. Although it has 
been shown (Goodman and LeVeque [6]) that TVD schemes are at most first order accurate 
in multiple space dimensions, it appears that entropy satisfying schemes do not share this 
restriction. Furthermore they are, in principle, extensible to the Navier-Stokes equations 
with finite rate chemistry. 

2. Generalized Entropy and Convection Processes 

Convection processes are often described in terms of control volumes. For a fixed 
control volume, a general convection process can be written 

( 2 . 1 ) J q it d"V 4 - j> f(q) • da = 0 

where the integrals are over the volume and over the surrounding surface respectively. 
The vector dn is the usual outward facing unit normal. For the case of gasdynamics, 
this covers conservation of mass, momentum and energy. Satisfaction of the first law of 
thermodynamics follows as a direct consequence. Unfortunately Eq. 2.1 says nothing 
about the second law of thermodynamics and admits unphysical solutions such as expansion 
shocks. 

For uniqueness, the second law of thermodynamics must also be satisfied. This can 
also be stated for a control volume. 


(2.2) J S{q), t dV+f F(q)-n>0 

In this equation S(q) is the generalized entropy, a scalar function of q, and F (q) is 
the so called generalized entropy flux. The quantity F(g) is properly a vector in the sense 
that it has a magnitude and a direction but not in the sense of q which contains several 
independent quantities. In one dimension, no generality is lost by representing F as a scalar 
quantity. Though not limited to the thermodynamic definitions, 5 and F must be chosen 
to satisfy the following relations. 


2 


(2.3) 


S i<jq < 0 


(2-4) — S t qf q 

The comma notation used here denotes differentiation. In the case where the indepen- 
dent variable q is a vector of unknowns, Eq. 2.3 is taken to mean that S tqq is a negative 
definite matrix. If q is a vector with more than two unknowns, the existence of 5 and F is 
not guaranteed. In any case, 5 may not be unique. For the case of the Euler equations it 
can be verified by substitution that the thermodynamic entropy satisfies these equations. 
Further discussion of generalized entropy can be found in a classic paper by Lax [7] . The 
numerical analog of Eqs. 2. 2-2.4 will be referred to as the entropy condition. 

There is one other requirement. The fluxes f(q) must remain bounded everywhere. 
For example, the Kutta condition is a way of satisfying this requirement near an airfoil 
trailing edge. The numerical analog of this requirement will be called the flux consistency 
condition. The flux consistency condition used in this work is 


(2 5) </(«,+ l) - /(»,+ *))(/(«,■+*) - /<«,» > o 

which requires that extrema in the flux occur at grid points. This is slightly stronger than 
the consistency condition suggested by Lax [7]. Lax suggests 


( 2 - 6 ) /(<7; + i) = /(?. j + i) = filj) ^ 9* = 9j + i 

which amounts to consistency in a Taylor series sense. The value of a flux consistency 
condition is twofold. The first is accuracy. Satisfaction of Eq. 2.6 in a stable scheme is 
enough to ensure convergence to the correct weak solution as the grid is refined in space 
and time. The second is entropy satisfaction. Satisfaction of Eq. 2.5 is needed to ensure 
that the entropy inequalities are satisfied between gridpoints. 

For various historical reasons, and because the notation is more compact, Eqs. 2.1 and 
2.2 are often replaced by an equivalent partial differential equation (PDE). For example Eq. 
2.1 is equivalent to the PDE 


(2-7) q,t + /(g),, = 0 

in one dimension. The remainder of this work will use the simpler PDE notation, though 
in principle the control volume notation could have been used throughout. 

This work deals with numerical satisfaction of the entropy condition and the flux 
consistency condition. These often ignored areas are important because they make the 
solution unique, and because they are the qualitative reason given for artificial dissipation 
or “smoothing.” 
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3. Fun with the Wave Equation 


The simplest convection model available is the scalar wave equation. 


(3.1) u, t + cu, x = 0 c > 0 

This linear equation has a simple analytic solution and does not support shock waves 
in the usual sense. In spite of this simplicity, it is not trivial to produce a good numerical 
algorithm for the wave equation. A good numerical algorithm should be conservative, second 
order accurate in space, and should satisfy a cell entropy for any CFL number sufficiently 
small. 

The simplest second order scheme is explicit Euler in time, central differenced in space. 

(3.2) xt 3 = Uj - ^(u i+ i - ti,-i) 



In this notation (which is only suitable for schemes with two time levels) superscripts 
indicate values at the new time level and subscripts represent values at the old time level. 
In the notation of Sweby this may be written 


(3.3) u J = Uj - vA xij .1 - A*(|i/Au J+ i) 
where we use the convention 

A + u, = Au J+ 1 = A'uj + i = u J+1 - Uj 

As written, the scheme is known to be unconditionally unstable. Traditionally one 
adds a second order dissipation term. 

(3.4) u J = u, - i/Au r | - A*(|i/Au J+ i) + A'(f j/a J + i Au J+ i) 

The variable a J + 1 is a smoothing coefficient which may or may not vary from point 
to point. Combining terms and using + i = 1 - + gives 


(3.5) u 3 = u, - vAu : .i - A‘(|i/^ + iAu j+ i) 

If i is a constant, the scheme that results is the standard central difference scheme 
with second order artificial viscosity as proposed by Lax [7]. The coefficient of artificial 
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viscosity is a in this case. If i = 0, the first order upwind scheme results. If ^ - + ± = 1, 
Eq. 3.2 is recovered. If + i =2, a first order downwind scheme results, suitable if c < 0. 

All constant values of <f> 3 +i , except = 0, result in schemes with spurious os- 

cillations, otherwise known as wiggles. A wiggle free, second order scheme will require 
nonconstant values of <t>j+±- One way to produce such a scheme is with the TVD ideas 
formalized by Harten and explained by Sweby. 

Expanding the last term of Eq. (3.5) gives 


(3.6) 


.J — 


- Uy - vA Uy_l - Au 3+i + |^ r lAu r l 


The standard practice at this point is to introduce the quantity r 3 = 
3.6 as 


Au 


- A 


Au . , i 

i+* 


and write Eq. 


(3.7) 


u ] = u j - i/Au r i - {^Au ri + Auy.l 


This allows collection of terms so that Harten ’s TVD conditions [8] can be applied. 


(3.8) 


U J ~ u, — V 


4> 

1 ( 1 : 


1+ilMr 1 -* ri) 


Au.i 

J 2 


The TVD proof says that the scheme will be TVD if we can insure that 


(3.9) 


0 < 1/(1 + ±$ y ] < 1 


A. 

where = (— y ^ — ^.i). With a little algebra the expression can be written (for v > 0) 
as 


(3.10) -2 < Qj < — - 

Inequalities (3.10) are plotted in Fig. 1. Everything in the shaded region is TVD. 
Generally, $ 3 takes on both positive and negative values in the neighborhood of zero. The 
object is to pick a maximum CFL number such that a reasonable range of can be used. 
Since $ 3 is independent of v the object is to pick rectangular regions from the shaded region. 
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Figure 1. Schemes in the shaded region are TVD. 



Figure 2. Schemes in the shaded region are TVD for v < 

For example, if u >> 1 then = -2. Unfortunately, this leads to an inconsistent 
scheme, u J+ i = uy, which is stable and TVD but gives the wrong answer. Consistent 



schemes (in the sense of Eq. 2.5) are possible if u < 1. If v = | the constraints are 


(3.H) 


-2 < 


<t> 


J+2 


- < 2 


where we have substituted in the definition for $ 3 . Generally, the coefficients <A J + j L and 
4> r i should be picked independently for algorithmic simplicity. Assuming that <£ J + ± can 
be chosen as a function of r : , inequalities 3.11 can be satisfied within the shaded region of 
Fig. 2. Negative values of 4> r i imply schemes with more dissipation than the first order 
upwind scheme and are excluded. 

Once again the task is to pick a suitable rectangular region from the shaded part of 
Fig. 2. Except for extrema it is easily shown that r ; as 1, with equality being achieved in 
the limit as Ax — » 0. Since all numerical analysis is based on this limiting case, we had 
better include it in our plans. A similar argument holds that except for discontinuities, 
ry-i as ry so that <pj+i 4>j.i ■ This being the case most of the action will occur on a line 
of slope one through the origin. The rectangular region with the largest amount of this line 
is the square defined by 


(3.12) 


0 <( 




*i-*) < 2 



Figure 3. Schemes in which <j>(r) is within the shaded region are TVD. 
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It is important to see that these inequalities must be satisfied at each point and that 
4>j+± w dl replace <p ]m L at the point j + 1. Thus, <t> r i can be replaced with <f> } + x in Eq. 
(3.12) and the resulting inequalities plotted in Fig. 3. The entire line = 0 is in- 

cluded in the shaded region. Notice that <f> = 0 represents the only choice for a constant 
coefficient TVD scheme of this type, and it is first order accurate. The idea here is to 
construct some function <p(r) that lies in the shaded region and results in a second or- 
der scheme. The function <f> = 1 yields the standard central difference scheme with no 
dissipation. The function <f> = r yields a three point, second order upwind scheme. Nei- 
ther of these is stable without dissipation and neither is TVD. Notice that both schemes 
pass through the point <j> = r = 1. Recall that r — > 1 as Ax — ► 0 for all schemes. For 
accuracy it is desirable that <£(1) = 1 so that the dissipation will vanish as the mesh is 
refined. All functions <f>(r) that pass through this point will be second order except at 
extrema and discontinuities. One function that lies between the two linear second order 
schemes and passes through the point <j> = r = 1 is shown in Fig. 3. Its formula is 


(3.13) 


*i+4 = 0, r J <0 


i 

^ + 5 r,- + r 


r, > 0 


This limiter is attributed to van Leer, though intended for use on Lax-Wendroff type TVD 
schemes. This scheme is programmable and vectorizable, forming successively r : from the 
definition, <f> J+ i from Eq. 3.13, and finally u 3 from Eq. 3.5. 


0 COMPUTED 
EXACT 




Figure 4. Results for a conservative TVD scheme that violates entropy. 


Figure 4 shows the performance of this simple TVD scheme. On the square wave initial 
data at a CFL number of | the solution is nearly perfect, as good as the mesh will permit. 
The reason is clear when the case of sine wave initial data is examined. Here too, the 
solution is a square wave. In fact, for an arbitrary initial pulse this scheme will eventually 
produce a convecting square wave. 

Intuitively a scheme that produces nothing but square waves seems to be a bad scheme. 
What makes it a bad scheme? Not conservation, for u is conserved exactly (the area under 
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the pulse is unchanged. Not wiggles, the scheme is TVD, as advertised. It is the failure of 
this scheme to satisfy the entropy inequality, important even in the absence of shock waves 
and discontinuities, that results in unphysical solutions. 

4. Entropy and the Wave Equation 

The definitions for entropy and the associated entropy flux are not unique. One such 
pair is S — -u 2 and F = -cu 2 . It can be verified by substitution that this satisfies Eqs 
2.3 and 2.4 as required. 

The actual entropy inequality Eq. (2.2) can only be evaluated if something is known 
about how S(u) varies over the control volume. In this work u is assumed to be piecewise 
constant so S(u) is constant over the control volume. Under this assumption u J+ ± is 
undefined, taking on any value between u } and u J + j. 

Using the explicit Euler time advance as before gives the difference scheme 


( 41 ) u 3 = Uj - - u } . l) 

and the corresponding discrete entropy statement 

(4.2) (ttO a -u; + i/(«; +4 -«^)<o 



Figure 5. Only schemes in the shaded region satisfy the entropy condition. 

Squaring Eq. 4.1 and substituting into Eq. 4.2 gives a discrete entropy statement that 
may be explicitly calculated. Remarkably, this factors, even without specifying u J + i. 
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Figure 6. Some TVD schemes satisfy an entropy constraint. 


(4.3) 



u 7 ) 2 - (u, -u r i) 2 + t/(u J+ i - u } . 


) 2 < 0 


Assume that v is positive. (The case in which v is negative is a trivial extension.) It is 
apparent that the third term will dominate for large enough v , leading to a kind of CFL 
restriction on satisfaction of the entropy inequality. The nature of the restriction becomes 
apparent when Eq. 4.3 is plotted as in Fig. 5. The shaded region indicates satisfaction of 
the entropy inequality. As in the previous section, the task at hand is to select a maximum 
value of v at which to operate. Whatever value is selected, the following inequalities apply. 


(4.4) 


_!< + 

~ XI j - U 1 ” V + 1 


So what does all this have to do with TVD schemes? If the desired difference scheme 
is given by Eq. 3.5, then appropriate definitions are, as verified by substitution, 


(4.5) + i = u : + ^ J+ i(u J+ i - Uj) 

and an equivalent formula for Uy.i. With these definitions, inequalities (4.4) simplify to 


(4.6) 


-1 < 


2-^-i 


< 


v - 1 

V + 1 
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These inequalities are plotted in Fig. 6. The smaller the CFL number, the larger the 
entropy satisfying region. The axes are the same as for Fig. 2, allowing comparison of the 
TVD region with the entropy satisfying region. In this case entropy satisfying schemes are 
a subset of TVD schemes. 

It is also worth noting that entropy satisfaction alone is insufficient; flux consistency 
(in the sense of Eq. 2.5) is also required. In Fig. 6 there is an entropy satisfying region to 
the right of <j> — 2 that does not satisfy the flux consistency requirement. Schemes operating 
in this region are generally unstable. 

As explained in the previous section , the case r 3 — 1 is the most important. This 
leads to the choice of a rectangular region in Fig. 6 described by 

(4.7) 0< 4> r x. <1-1/ 

-1 - 1 / < < 1 - 1 / 

r. 



Figure 7. Schemes in the shaded region satisfy an entropy constraint. 


This may be compared with Eq. (3.12). Finally these inequalities are plotted in Fig. 
7. The shaded region represents a sufficient condition for entropy satisfaction in this case. 
For comparison, the TVD region is also shown, as is the limiter derived in the previous 
section. Notice that it lies outside the entropy satisfying region. 

On the other hand, it is possible to conclude that the limiter represented by the edge 
of the entropy satisfying region, will definitely lead to a scheme that satisfies the entropy 
condition on a cell by cell basis. From the previous section it is easily shown to be a TVD 
scheme as well. Unfortunately it does not, in general, contain the point <f> = r 3 = 1 as 
required for second order accuracy unless v = 0. It is possible to show (by analysis not 
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included here ) that such a scheme has spatial accuracy of order At -+- Az 2 . This means 
that it is a first order scheme unless At is of the order of Az 2 . The spatial accuracy, and 
hence the steady state solution, depends on the time step, a feature that this approach 
was supposed to avoid. Numerical results are shown in Fig. 8. As expected, the scheme 
is highly dissipative for a CFL number of It is, however, less dissipative than the first 
order upwind scheme (and a lot more complicated). That’s progress. 


0 COMPUTED 
EXACT 
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Figure 8. Results for a first order, entropy satisfying scheme. 


5. Implicit schemes 

Any time a scheme, such as the one just described, shows good features at small CFL 
numbers that are lost at large ones, it may be possible to avoid the problem using an implicit 
scheme. In this case an implicit Euler scheme was selected to make use of all the analysis 
in Sections 3 and 4. 


(5.1) u J = u } - ~(u J+1 - u 1 ' 1 ) 

£ 

Collecting terms according to time level gives 

(5.2) u J + _ u Jm *) — Uj 

Recall that the entropy statement must be differenced (integrated) the same way as the 
PDE. This leads to a slightly different entropy statement: 

(5.3) (u J ) 2 - u 2 + t/((«J + = ) 2 - (u J **) 2 ) < 0 
Squaring Eq. 5.2 and substituting for u* in Eq. 5.3 gives, for positive v, 
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Figure 9. Entropy satisfying regions for implicit scheme. 


(5.4) (ti J+ 5 - u 3 ) 2 — (u J - u 3 ~ - ) 2 — 2 — u 3 ~ 2 ) 2 < 0 

Notice the similarity to Eq. 4.3 . Except for the time level, the only difference is in 
the sign of the third term. To see what effect this has, inequalities 5.4 are plotted in Fig. 9. 
The shaded region represents satisfaction of both inequalities. This region is much larger 
than the equivalent region for explicit Euler time advance (Fig. 5). The entropy inequality 
will be satisfied, independent of v , so long as 


(5.5) 


-1 < 


U J+2 - U 3 
U 3 — U 3 ~ 2 


< 1 


Following the methodology shown in the previous section leads to the limiter 


(5.6) <f> 3 + 2 =mm(l,|r 7 |) 

which is plotted in Fig. 10. Except for the difference in time level, this is the same limiter 
computed in the previous section if v — 0. Its advantage occurs for larger values of v. 
Notice that central differencing is used except for a small region around the origin. This 
means that the accuracy of the scheme may (or may not) remain second order accurate 
at extrema (negative values of r 3 ). By contrast, the limiter shown in Fig. 3 will always 
degenerate to first order at extrema. 

There is still the problem of computing tz J+ 2 at the new time step. This can be approx- 

1 du + x , 

imated by a Taylor series type linearization. For example + x -f — (u* “ u k ) 
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Figure 10. Limiter for implicit TVD scheme. 
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Figure 11. Results for implicit scheme. 


where u J+ i is computed just like u J+ 2 , but at the old time level. The k subscript is 
summed, as in tensor notation, over all values of j (most of the terms vanish). A system of 
linear equations can be formed that has two bands below the diagonal and one above the 
diagonal. In a more general situation this would be a pentadiagonal. In computing the par- 
tials, the only difficulty comes in linearizing the limiter, <p(r). This must be done piecewise 
and can be a poor approximation in the transient at high CFL numbers. Numerical results 
are plotted in Fig. 11. At low CFL numbers the scheme performs very well, producing 
a TVD result that satisfies the entropy inequality. At high CFL numbers [v >> 1) the 
entropy inequality (and the TVD condition) may not be satisfied in the transient because 
of errors introduced by the linearization. The scheme remains stable but becomes very 
diffusive due to the first order accurate time advance. At CFL numbers between u — .5 and 
v = 1 the scheme may actually perform worse than the first order upwind scheme. This is 
due to the latter schemes special properties for the wave equation at v = 1. 
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6. Nonlinear Schemes 


It is reasonable to suggest that the wave equation is irrelevant, that it is too far 
removed from the interesting applications. In spite of this, few schemes have good analysis 
for anything else. In this case, however, more complete analysis is possible. The case of 
scalar, nonlinear equations is covered in this section. 

Consider the scalar, nonlinear , convection equation. 


(61) ?,t + /(?), x = 0 

For this section, at least, q and / are considered to be scalars. For example, the special 
case f(q) = ^q 2 is often known as Burger’s equation. As in the last section, an implicit 
difference scheme is selected. 

(6.2) ^- 9j + A(/;+5-/H) = 0 

Instead of u the definition A = is used since the wave speed is no longer constant. 

The entropy statement in this case is 


(6.3) S 3 - Sj + A (F 3 + i - F j -*)> 0 

Since the generalized entropy S can be a somewhat complicated function of q, it may not be 
practical to simply substitute Eq. 6.2 into 6.3 as was done for the wave equation. Instead, 
Sj is expanded in a truncated Taylor series. 

(6.4) S } = S 3 -S 3 q Aq+\S, qq ^J)Aq 2 


In Eq. 6.4, Aq is defined by Aq = q 3 - qj. Also, the third term is evaluated at a time 
£, an unspecified time between the old and new time levels whose existence is guaranteed 
by the intermediate value theorem. Equation 6.4 is exact, all the uncertainty is tied up in 
the value of f. Substituting this relation into Eq. 6.3 gives 


(6.5) S 3 q Aq^X(F 3 ^- - F 3 ^) > ±S, qq (Z,j)Aq 2 


Since S(q) must satisfy Eq. 2.3 (convexity), the right hand side of Eq. 6.5 is always 
negative or zero. It is zero at steady state (Ag = 0). A more restrictive condition (which is 
easier to implement than Eq. 6.5) is 
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( 6 . 6 ) 


S’ q Aq + \{F ] + ± - F J '*) > 0 

Notice that Eq. 6.2 gives a relation for A q. Substituting this into Eq. 6.6 and dividing out 
A gives 


(6.7) -5i(/ J ' + * - f 3 '*) + (F’ + * - F r *) > 0 

Equation 6.7, though correct, is difficult to implement because of terms evaluated 
between grid points. The solution selected here is to expand such terms in a Taylor series 
about the point j to two terms. For example the approximation 


( 6 . 8 ) r^Kf’ + f^&q+y^Aq 2 

is used. After dropping the obvious terms, what remains is 


(6.9) («+ t)(F J , - Sj/’,) + i(a ! - - Sj/j,) > 0 

where a = q 3+ ? - q } and b = q } - q 3 ~% . The first term is zero, since S(q) and F(q) must 
satisfy Eq. 2.4. The second term simplifies substantially using 


(6- 10) F qq — (F q ) q — (S q f q ) q — S <qq f iq + S iq f iqq 

which comes from using Eq. 2.4 and the chain rule. Using Eq. 6.10 to simplify Eq. 6.9 
gives the entropy statement, 


(6-11) l 2 (a 2 -b 2 )S tqq f, q >0 

which simplifies further using 5 q<? < 0. The final result is 


(6.12) sign(/ i9 )(a 2 /6 2 - 1) < 0 

The quantity f <q is just the local wave speed. The term sign (f q ) appears to be re- 
quired for satisfaction of the discrete entropy inequality. Instead of being postulated for 
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physical reasons (as in upwind schemes), it appears directly in the algebra. The linearization 
difficulties posed by such a term are considerable. 

In spite of the simplicity of Eq. 6.12, there are a few complications to solving nonlinear 
problems. First, the discrete entropy inequality is not satisfied exactly, as it was for the 
wave equation, but only to second order. This is a consequence of approximations such as 
Eq. 6.8. Second, and more serious, is the problem of sign changes in the local wave speed 
(shocks and sonic points). 

Equation 6.12 must be satisfied at each node point. The value of q 3 + 2 that is chosen 
must be compatible with inequalities at point j and at point j + 1. The important feature 
to take into account is the signs of the wave speeds at these two points. 

The case in which both are positive was covered in the previous section. Equation 6.12 
is exact for the wave equation and leads to exactly the same algorithm previously discussed. 
The case in which both are negative is easily derived by reversing the coordinate system. 
If, by chance or design, / i9 = 0 for a particular node, Eq. 6.12 will always be satisfied. The 
interesting cases are the sign changes, the shock, and sonic point cases. 

A sonic point is said to occur if f 3 q < 0 and f 3 q +1 > 0. The characteristics are said to 
be diverging at this point. This is a problem since 6 : depends in general on the ratio r, 
or r J + i depending on which is the upstream point. At the sonic point neither side can be 
said to be upstream. In choosing <j> 3 '^ it is subtly assumed that 1 < <t> 3 + - <2. In choosing 
4> 3+ i it is subtly assumed that 0 < <p 3 + 2 < 1. The only solution that satisfies entropy 
inequalities for cells j and j + 1 is <)> 3+ * = 1. 

A shock point is said to occur if f 3 q > 0 and f 3 + l < 0. The characteristics are said to 
be converging at this point. This case is the most difficult. Since q(x) is discontinuous at 
the shock, Taylor series approximations such as Eq. 6.8 are not valid. Equation 6.12 must 
be replaced by two equations, one on each side of the shock 

(6.13) (F 3+ *- - F 3 ) - S 3 q (f 3+ '- - f 3 ) > 0 

( F 3+1 - F 3+ =) - S 3 q +1 {f 3+1 - f 3 + '-) > 0 

The endpoint values, <j> — 0 and <f> = 2, make good guesses because they at least satisfy 
one of equations 6.13. For the special case of Burger’s equation it can be shown that they 
are the only possible solutions. Furthermore, for some combinations of q 3 and + 1 only 

one of these is valid, 


(6.14) <^ + = = 1 - sign(/ J ? + f 3 q +1 ) 

This says that on an interval with a shock in it, 2 — u J or + The choice 

is determined by the sign of the shock speed, which is consistent with characteristic theory. 
Assuming that Eq. 6.14 is valid in general, the limiter formula’s for a scalar nonlinear 
equation are summarized in Table 1. 
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Table 1 - The Function <f>(r) 


<f>} + 2 = 

P, q < 0 

f% > o 

f 3 q +1 < o 

2 - min(l, |r J + 1 |* 1 ) 

1 - sign(/ J ? + //,+ 1 ) 

/^ +1 > 0 

1 

min(l, |r J |) 


For the special case of Burger’s equation, a suitable choice of entropy and entropy flux 
is S(q) = q 2 , F(q) — | q 3 . An interesting initial condition is q(x) — sin(i) on the interval 
x — [0, 2tt] - This case was computed for several values of A. Figure 12 shows the results 
for an initial CFL value of The computed values are shown as symbols whereas the 
exact values are shown as a solid line. Also shown is <f>( x) at the beginning and end of 
the calculation. Notice that central differencing is used in expansion regions and second 
order upwind differencing is used in compression regions. At the shock a first order upwind 
scheme is employed. Once the shock has formed, all of the dissipation occurs at the shock, 
a feature of the PDE. 



Figure 12. The differencing is a function of the solution. 


There are a few anomalies, not apparent from Fig. 12. The limiters are computed 
in such a way that only a single interval is allowed for the shock. Yet, in this worst case 
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scenario, the shock lands on a grid point. What appears to be a two interval shock is 
treated by the scheme as a shock and a strong compression. If, due to roundoff, the shock 
is initially moving to the right, the center point will rise above zero. This will result in 
the shock moving to the right on the next step. After the solution matures, there is no 
visible difference between the shock and the compression. Another problem is that pure 
central differencing is neutrally stable. Occasionally, rounding errors will result in a visible 
excursion of the limiter away from central differencing even in the smooth regions. These 
excursions are very small, however, and the solution is not noticeably affected. 

At moderate CFL numbers (A « 1) the scheme suffers from minor time accuracy prob- 
lems in the transient. This is a result of the first order time advance. At high CFL numbers 
(A >> 1), the scheme loses the TVD property in the transient because of linearization 
errors. It remains stable, however, and achieves the correct steady state solution at CFL 
numbers of 1000. 


7. Conclusions 

Initially this work was undertaken to explain artificial dissipation (“smoothing”) in 
physical terms. The explanation is that smoothing is a numerical artifice to satisfy the 
discrete entropy inequality in a global sense. This ensures stability (given suitable boundary 
conditions) since the total entropy within the domain must be non-diminishing. Entropy 
in a region reaches a maximum when the solution is a constant over the entire domain. 
Normally a balance is reached before this happens, but in any case the solution cannot 
blow up. If a global entropy inequality is satisfied while a local one is not, there will be 
several cells (for example, near shocks) that destroy entropy. This allows the oscillations 
so familiar with linear smoothing. Physically, the solution is only meaningful if integrated 
over a collection of cells that, together, satisfy the entropy inequality. Near shocks, this 
may require several cells to be lumped together and averaged. Thus the effective grid near 
a shock may be coarser than it appears if the smoothing does not satisfy a cell entropy 
inequality. 

It should also be mentioned that smoothness of a solution is no guarantee that a local 
entropy inequality is satisfied. The author has observed smooth solutions of the quasi-one- 
dimensional Euler equations for which every cell in the supersonic region consumes entropy. 
Such solutions do not correspond to flows of very viscous fluid as is often thought. That 
they do not correspond to any physical fluid is evidenced by the fact that the second law 
of thermodynamics is violated over much of the domain. 

The flux consistency condition puts bounds on what happens between grid points. In 
this work a piecewise constant assumption is used to compute the discrete entropy inequality. 
The evaluation of flux between gridpoints must be consistent with this assumption. The 
flux consistency condition used here is not unique, just sufficient. It is only valid if the 
piecewise constant assumption is used. 

The purpose of this work is not so much to recommend a particular scheme as it is to 
recommend a method of analysis based on satisfying the entropy inequality and a certain 
flux consistency condition. The scheme derived here is offered as proof that schemes which 
satisfy these two conditions do exist. 
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This work has shown that the TVD property is insufficient to guarantee a physical 
solution, either in the steady state or in the transient. To provide such a guarantee, a 
scheme must satisfy the discrete entropy inequality for time dependent problems, or the 
semi-discrete entropy inequality for steady state problems. Schemes that satisfy a discrete 
entropy inequality may need to be implicit to achieve second order accuracy in space. 
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